Recap: we have one sample

library(patchwork)
library(ggplot2)

set.seed(108)
heights <- rnorm(100, mean = 170, sd = 10)
heart <- rnorm(100, mean = 70, sd = 10)

# heights
p1 <- ggplot(data.frame(heights), aes(x = heights)) +
  geom_histogram(aes(y = after_stat(density)),
    binwidth = 5, fill = "lightblue", color = "black") +
  geom_vline(aes(xintercept = mean(heights)), 
    color = "red", linetype = "dashed") +
  ggtitle("Height of students") +
  theme_classic()


# heart rate
p2 <- ggplot(data.frame(heart), aes(x = heart)) +
  geom_histogram(aes(y = after_stat(density)), 
    binwidth = 5, fill = "lightblue", color = "black") +
  geom_vline(aes(xintercept = mean(heart)), 
    color = "red", linetype = "dashed") +
  ggtitle("Heart rate") +
  theme_classic()


p1 + p2

What if we want to compare a sample of data to another sample?

set.seed(129)
## Generate data for 2 groups
heights_all <- data.frame(
  group = rep(c("ENVX1002", "ENVX2001"), each = 100),
  heights = c(rnorm(100, mean = 165, sd = 7), 
    rnorm(100, mean = 185, sd = 9)))
heart_all <- data.frame(
  group = rep(c("ENVX1002", "ENVX2001"), each = 100),
  heart_rates = c(rnorm(100, mean = 70, sd = 10), 
    rnorm(100, mean = 72, sd = 5)))

## Plot
p3 <- ggplot(heights_all, aes(x = heights, fill = group)) +
  geom_histogram(aes(y = ..density..), 
    binwidth = 5, color = "black", alpha = .2) +
  geom_vline(aes(xintercept = 170), 
    color = "red", linetype = "dashed") +
  geom_vline(aes(xintercept = 185),
    color = "blue", linetype = "dashed") +
  ggtitle("Height of students (cm)") +
  theme_classic()
p4 <- ggplot(heart_all, aes(x = heart_rates, fill = group)) +
  geom_histogram(aes(y = ..density..), 
    binwidth = 5, color = "black", alpha = .2) +
  geom_vline(aes(xintercept = 70), 
    color = "red", linetype = "dashed") +
  geom_vline(aes(xintercept = 72), 
    color = "blue", linetype = "dashed") +
  ggtitle("Heart rates (bpm)") +
  theme_classic()

p3 + p4

Comparing two samples: visualisation

p5 <- ggplot(heights_all, aes(x = group, y = heights, fill = group)) +
  geom_boxplot(alpha = .2) +
  ggtitle("Height of students (cm)") +
  theme_classic()
p6 <- ggplot(heart_all, aes(x = group, y = heart_rates, fill = group)) +
  geom_boxplot(alpha = .2) +
  ggtitle("Heart rates (bpm)") +
  theme_classic()

p5 + p6

Some considerations: the boxplot

Data

A simulated example (data is not real):

redbull <- data.frame(
  group = c(rep("redbull", 12), rep("control", 12)), 
  heart_rate = c(72, 88, 72, 88, 76, 75, 84, 80, 60, 96, 80,  84, 84, 76, 68, 80, 64, 62, 74, 84, 68, 96, 80, 64))

str(redbull)

'data.frame':   24 obs. of  2 variables:
 $ group     : chr  "redbull" "redbull" "redbull" "redbull" ...
 $ heart_rate: num  72 88 72 88 76 75 84 80 60 96 ...

Hypothesis

For a two-sample t-test, the null hypothesis is that the means of the two groups are equal, and the alternative hypothesis is that the means are different.

H_0: \mu_{\text{redbull}} = \mu_{\text{control}} H_1: \mu_{\text{redbull}} \neq \mu_{\text{control}}

Assumptions

The assumptions of the two-sample t-test include:

1. Normality: the data are normally distributed.
2. Homogeneity of variance: the variances of the two groups are equal.

Normality: histogram

• We visually inspect the distribution of the data using histograms, generally for each group.
• Look out for: symmetry, skewness, and multimodality.
• Hard to visualise when n (sample size is small)

par(mfrow = c(1, 2))
hist(redbull$heart_rate[redbull$group == "redbull"], 
  main = "Red Bull group", xlab = "Heart rate (bpm)", 
  col = "lightblue", border = "black")
hist(redbull$heart_rate[redbull$group == "control"],
  main = "Control group", xlab = "Heart rate (bpm)", 
  col = "lightblue", border = "black")

library(ggplot2)
ggplot(redbull, aes(x = heart_rate, fill = group)) +
    geom_histogram(aes(y = ..density..),
        binwidth = 12, color = "black", alpha = .5
    ) +
    facet_wrap(~group) +
    theme_classic()

library(ggpubr)
gghistogram(redbull, x = "heart_rate", fill = "group", binwidth = 6) +
    facet_wrap(~group) +
    theme_classic()

Normality: QQ-plot

• The qq-plot is a graphical method to specifically assess the normality of the data. Again, we look at the data for each group.
• Look out for: deviations from the straight line.

par(mfrow = c(1, 2))
qqnorm(redbull$heart_rate[redbull$group == "redbull"],
    main = "Red Bull group", xlab = "Theoretical quantiles",
    ylab = "Sample quantiles", col = "blue"
)
qqline(redbull$heart_rate[redbull$group == "redbull"], col = "red")

qqnorm(redbull$heart_rate[redbull$group == "control"],
    main = "Control group", xlab = "Theoretical quantiles",
    ylab = "Sample quantiles", col = "blue"
)
qqline(redbull$heart_rate[redbull$group == "control"], col = "red")

library(ggplot2)
ggplot(redbull, aes(sample = heart_rate)) +
    stat_qq() +
    stat_qq_line() +
    facet_wrap(~group) +
    theme_classic()

library(ggpubr)
ggqqplot(redbull$heart_rate) +
    facet_wrap(~ redbull$group) +
    theme_classic()

Normality: formal test

• Use the Shapiro-Wilk test which tests the null hypothesis that the data are normally distributed.
• This test is sensitive to deviations from normality in the tails of the distribution, and is suitable for small sample sizes (about 5 to 50 observations).

shapiro.test(redbull$heart_rate[redbull$group == "redbull"])

    Shapiro-Wilk normality test

data:  redbull$heart_rate[redbull$group == "redbull"]
W = 0.97459, p-value = 0.9524

shapiro.test(redbull$heart_rate[redbull$group == "control"])

    Shapiro-Wilk normality test

data:  redbull$heart_rate[redbull$group == "control"]
W = 0.93733, p-value = 0.4643

What if the normality assumption is violated?

• The t-test is robust to deviations from normality, especially for large sample sizes due to the Central Limit Theorem.
• If the sample size is small, consider using a non-parametric test (e.g. the Wilcoxon rank-sum test): next week
• Alternatively, transform the data: later

Equal variances: boxplot

• We visually inspect the spread of the data using boxplots, generally for each group.
• Look out for: differences in spread, outliers, and symmetry.

par(mfrow = c(1, 2))
boxplot(heart_rate ~ group,
    data = redbull,
    main = "Heart rate", xlab = "Group", ylab = "Heart rate (bpm)"
)

library(ggplot2)
ggplot(redbull, aes(x = group, y = heart_rate, fill = group)) +
    geom_boxplot(alpha = .2) +
    labs(x = "Group", y = "Heart rate (bpm)") +
    theme_classic()

library(ggpubr)
ggboxplot(redbull,
    x = "group", y = "heart_rate",
    fill = "group", alpha = .2
) +
  labs(x = "Group", y = "Heart rate (bpm)") +
  theme_classic()

Equal variances: formal tests

• Bartlett’s and Levene’s tests may be used to test the null hypothesis that the variances of the groups are equal.
• These tests are sensitive to deviations from normality (Levene’s test is less so compared to Bartlett’s), and are suitable for small sample sizes.

library(car)
leveneTest(heart_rate ~ group, data = redbull)

Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  1   0.289 0.5962
      22               

bartlett.test(heart_rate ~ group, data = redbull)

    Bartlett test of homogeneity of variances

data:  heart_rate by group
Bartlett's K-squared = 0.075121, df = 1, p-value = 0.784

What if the equal variance assumption is violated?

Some debate exists on what to do, but choices include:

• Use the Welch’s t-test, which is robust to unequal variances: coming up next
• Transform the data to stabilise the variances: later
• Perform a non-parametric test: next week

The Welch’s t-test

• The Welch’s t-test is a modification of the two-sample t-test that does not assume equal variances.
• Also applicable when the sample sizes are unequal.

Why not use the Welch’s t-test all the time?

• Ongoing debate on whether to use the Welch’s t-test or the Student’s t-test when the variances are equal.
• The Welch’s t-test is generally considered more robust, and is the default in R’s t.test() function.
• You can still use the Student’s t-test by setting var.equal = TRUE in the t.test() function.

Are the assumptions of normality and homogeneity of variance met?

When reporting in journals, it is common to simply state that the assumptions were met and what tests were used to confirm them, without showing the exact results of the tests!

Performing the t-test

t.test(heart_rate ~ group, data = redbull, var.equal = TRUE)

    Two Sample t-test

data:  heart_rate by group
t = -1.1365, df = 22, p-value = 0.268
alternative hypothesis: true difference in means between group control and group redbull is not equal to 0
95 percent confidence interval:
 -12.947117   3.780451
sample estimates:
mean in group control mean in group redbull 
             75.00000              79.58333 

Results indicate that the means of the two groups are not significantly different (p = 0.27).

t.test(heart_rate ~ group, data = redbull)

    Welch Two Sample t-test

data:  heart_rate by group
t = -1.1365, df = 21.845, p-value = 0.2681
alternative hypothesis: true difference in means between group control and group redbull is not equal to 0
95 percent confidence interval:
 -12.950568   3.783902
sample estimates:
mean in group control mean in group redbull 
             75.00000              79.58333 

Conclusion

Differences in heart rate we not statistically significant between the Red Bull and control groups (t₂₂ = -1.1, p = 0.27) indicating that Red Bull did not significantly increase the heart rate of students sampled from ENVX1002.

Are the two sample independent?

When testing if two samples are different from each other, we need to consider two possible scenarios:

Paired t-test

Experimental design (what if?)

Before

Two groups of students selected at random, without replacement, from the ENVX1002 cohort.

Hypothesis

For a paired t-test, the null hypothesis is that the mean difference between the two groups is zero, and the alternative hypothesis is that the mean difference is different from zero.

H_0: \mu_{\text{diff}} = 0 H_1: \mu_{\text{diff}} \neq 0

Assumptions of the paired t-test

• The assumption of the paired t-test is that the differences between the two groups are normally distributed.
• There is no assumption of equal variances, as the paired t-test is a one-sample t-test on the differences.
  • Another way to think about it is that since the data are paired, the variance of the differences is the same for both groups.

Performing the paired t-test

There are two ways.

diff <- redbull$heart_rate[redbull$group == "redbull"] - 
  redbull$heart_rate[redbull$group == "control"]
t.test(diff)

    One Sample t-test

data:  diff
t = 1.6578, df = 11, p-value = 0.1256
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 -1.501628 10.668294
sample estimates:
mean of x 
 4.583333 

# t.test(heart_rate ~ group, data = redbull, 
#   paired = TRUE)

Recap: Assumptions of the two-sample t-test

With independent samples:

• Normality: the data are normally distributed
• Homogeneity of variance (equal variances): the variances of the two groups are equal

With paired samples:

• Normality: the differences between the paired samples are normally distributed
• Equal variances is implied

If we analyse the data anyway…

The t-test:

• may provide incorrect results as mean and variance calculations depend on normally distributed data.
• may be less powerful (i.e., less likely to detect a true difference).
• may be biased (i.e., systematically over- or under-estimating the true difference).

What can we do?

The t-test is quite robust to violations of normality, especially when the sample size is large. However, the assumption of equal variances is more critical – we cannot simply depend on large sample sizes to “fix” the problem.

Options include:

• Transform the data to normalise the data and/or scale the variance
• Use a Welch’s t-test or a Welch’s ANOVA (limited cases)
• Use a non-parametric test, such as the Mann-Whitney U test or Wilcoxon signed-rank test (paired samples) – however, these tests have less power than the t-test i.e. less likely to detect a true difference.

Is the food collected by ants different between two sites?

Data structure

library(tidyverse)
ants <- read.csv("data/ants.csv") %>%
  mutate(Tree = factor(Tree))

glimpse(ants)

Rows: 54
Columns: 2
$ Food <dbl> 11.9, 33.3, 4.6, 5.5, 6.2, 11.0, 24.3, 20.7, 5.7, 12.6, 10.2, 4.7…
$ Tree <fct> Rowan, Rowan, Rowan, Rowan, Rowan, Rowan, Rowan, Rowan, Rowan, Ro…

We want to compare the mean biomass of food, collected by ants between the two sites in dry weight (mg) of prey, divided by the total number of ants leaving the tree in 30 minutes.

Visualising the data

library(ggplot2)
p_ants <-
  ggplot(ants, aes(x = Tree, y = Food)) +
  geom_boxplot() +
  ylab("Biomass of food (mg per ant)") +
  theme_minimal()

p_ants

Does this data meet the assumptions of the two-sample t-test?

Checking assumptions

We have some idea that the data may not be normally distributed, but are not quite sure. So let’s check using the Q-Q plot.

ggplot(ants, aes(sample = Food)) +
  stat_qq() +
  stat_qq_line() +
  facet_wrap(~Tree) +
  theme_minimal()

Picking a transformation

We need to consider the type of data and the shape of its distribution when choosing a transformation. These can be assessed using:

Skewness

The degree of asymmetry in the data distribution when compared to a normal distribution.

Example: skewness

library(tidyverse)
library(patchwork)

x <- seq(0, 1, length.out = 100)

# Calculate the density of the Beta distribution at these points
data1 <- data.frame(x = x, y = dbeta(x, 5, 2), dist = "Negative (left) skewed")
data2 <- data.frame(x = x, y = dbeta(x, 5, 5), dist = "Symmetrical")
data3 <- data.frame(x = x, y = dbeta(x, 2, 5), dist = "Positive (right) skewed")
data <- rbind(data1, data2, data3) %>%
  mutate(dist = factor(dist, levels = c("Negative (left) skewed", "Symmetrical", "Positive (right) skewed")))

# Create the plot
ggplot(data, aes(x = x, y = y)) +
  geom_line(color = "skyblue") +
  geom_area(fill = "skyblue", alpha = 0.4) +
  theme_minimal() +
  theme(axis.text = element_blank(), axis.ticks = element_blank()) +
  facet_wrap(~dist) +
  ylab("density") +
  xlab("")

Kurtosis

Used to describe the extreme values (outliers) in the distribution versus the tails.

Example: kurtosis

library(tidyverse)
library(patchwork)
library(moments)

# Generate data
set.seed(123)
x1 <- seq(-10, 10, length.out = 1000)
x2 <- seq(-5, 5, length.out = 1000)
data1 <- data.frame(x = x1, y = dt(x1, df = 1), dist = "High Kurtosis")
data2 <- data.frame(x = x2, y = dt(x2, df = 10), dist = "Low Kurtosis")
data <- rbind(data1, data2) %>%
  mutate(dist = factor(dist, levels = c("High Kurtosis", "Low Kurtosis")))

# Create the plot
ggplot(data, aes(x = x, y = y)) +
  geom_line(color = "skyblue") +
  geom_area(fill = "skyblue", alpha = 0.4) +
  theme_minimal() +
  theme(axis.text = element_blank(), axis.ticks = element_blank()) +
  facet_wrap(~dist, scales = "free_x") +
  ylab("density") +
  xlab("")

Skewness and kurtosis in the ants data

With experience we can “eyeball” the data, but we can also calculate the skewness and kurtosis.

ants %>%
  group_by(Tree) %>%
  summarise(skewness = skewness(Food), kurtosis = kurtosis(Food))

# A tibble: 2 × 3
  Tree     skewness kurtosis
  <fct>       <dbl>    <dbl>
1 Rowan       1.04      3.15
2 Sycamore    0.807     2.63

Workflow

Picking a transformation

How do we check if the transformation worked?

We need to apply the transformation to the entire dataset and check the Q-Q plot again.

ants$Food_log <- log(ants$Food)

# compare the Q-Q plot before and after transformation
pfood <- ggplot(ants, aes(sample = Food)) +
  stat_qq() +
  stat_qq_line() +
  facet_wrap(~Tree) +
  ggtitle("Before transformation") +
  theme_classic()

pfoodlog <- ggplot(ants, aes(sample = Food_log)) +
  stat_qq() +
  stat_qq_line() +
  facet_wrap(~Tree) +
  ggtitle("After transformation") +
  theme_classic()

pfood / pfoodlog

Checking skewness and kurtosis after transformation

ants %>%
  group_by(Tree) %>%
  summarise(skewness = skewness(Food_log), kurtosis = kurtosis(Food_log))

# A tibble: 2 × 3
  Tree      skewness kurtosis
  <fct>        <dbl>    <dbl>
1 Rowan    0.271         1.77
2 Sycamore 0.0000457     1.86

Performing the t-test

fit <- t.test(Food_log ~ Tree, data = ants, var.equal = TRUE)
fit

    Two Sample t-test

data:  Food_log by Tree
t = -2.0521, df = 52, p-value = 0.04521
alternative hypothesis: true difference in means between group Rowan and group Sycamore is not equal to 0
95 percent confidence interval:
 -0.732858080 -0.008203447
sample estimates:
   mean in group Rowan mean in group Sycamore 
              2.287756               2.658287 

Back-transforming the results

• For power transformations, we can back-transform the results to the original scale using the inverse function.
• Log transformations are a bit tricky as the inverse function is the exponential function.
  • For the natural log transformation which is log() in R, the inverse function is the exponential function: e^x.
  • For the base 10 log transformation which is log10() in R, the inverse function is 10^x.

Interpretation

Back-transforming mean values

browan <- exp(fit$estimate[[1]]) # mean biomass from the Rowan site
bsycamore <- exp(fit$estimate[[2]]) # mean biomass from the Sycamore site

# check the ratio
bsycamore / browan

[1] 1.448503

Evidence suggests that the log-transformed mean biomass of food collected by ants from the Rowan site is significantly different from the log-transformed mean biomass of food collected by ants from the Sycamore site (t = -2.05, df = 52, p = 0.045).

The mean biomass of food collected by ants from the Sycamore site (14.3 mg) is 1.4 times greater than the mean biomass of food collected by ants from the Rowan site (9.9 mg).

Back-transforming confidence intervals

ant_ci <- exp(fit$conf.int)
ant_ci

[1] 0.4805336 0.9918301
attr(,"conf.level")
[1] 0.95

Comparing to a test without transformation

fit2 <- t.test(Food ~ Tree, data = ants, var.equal = TRUE)
fit2

    Two Sample t-test

data:  Food by Tree
t = -1.9217, df = 52, p-value = 0.06013
alternative hypothesis: true difference in means between group Rowan and group Sycamore is not equal to 0
95 percent confidence interval:
 -10.4916030   0.2267678
sample estimates:
   mean in group Rowan mean in group Sycamore 
              12.27143               17.40385 

References

• Quinn G. P. & Keough M. J. (2002) Experimental design and data analysis for biologists. Cambridge University Press, Cambridge, UK.
• Logan, M. (2010). Biostatistical design and analysis using R a practical guide. Hoboken, N.J., Wiley-Blackwell.